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ABSTRACT 


The Maximum Likelihood estimators for the Neyman Type A 
distribution parameters are very difficult to compute. [In 
this thesis, the Empirical Probability Generating Function 
is used to provide estimators that are easier to compute and 
have asymptotic efficiency at least as high as 97% of that 
for the Maximum Likelihood estimators over most of the param- 
eter space considered. The estimators found by this method 
are consistently better than the Method of Moments and the 
Method of Zero Frequency estimators with respect to asymptotic 
efficiency. The considerations of preference in using one 


method over another are discussed. 
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eee NU ne Duel LON 


Several 'contagious' distributions, including the Neyman 
Type A distribution, were derived by Neyman [1] in the course 
of modeling the reproduction behavior of certain types of 
plants and bacteria. More currently, these models come under 
the heading of branching processes. Several methods of pa- 
rameter estimation have been proposed for the Neyman Type A 
distribution with varying degrees of success. The two most 
frequently used are the Method of Moments and the Method of 
Zero Frequency. The performance of an estimation scheme can 
be measured by the asymptotic efficiency of its estimators 
as defined by Wilks [2]. This paper proposes a method for 
estimating the parameters of the Neyman Type A distribution 
by uSing the Empirical Probability Generating Function. 

For the Neyman Type A distribution, the Method of Moments 
provides estimates that are quick and easy to compute, al- 
though the efficiency of these estimators is not always high. 
It is well known that Maximum Likelihood estimators are 
asymptotically efficient, but they are very difficult to com- 
pute for the Neyman Type A distribution. 

The method of estimation by using the Empirical Probabil- 
ity Generating Function provides estimators with asymptotic 
efficiency close to that of the Maximum Likelihood estimators 


over most of the parameter space considered. Extensive tables 


a 





are presented along with an algorithm that is relatively easy 
to use on current programmable calculators. 

The estimators found by applying this method have con- 
Sistently higher asymptotic efficiency than the Method of 
Moments and the Method of Zero Frequency estimators. However, 
larger sample sizes may be required to cause the algorithm to 
converge on a unique pair of estimates. 

A brief outline of the report is presented as follows. A 
derivation of the Neyman Type A distribution with a list of 
some of its mathematical properties is contained in section 
II. The equations for the Maximum Likelihood estimators and 
asymptotic efficiency appear in section III. The mathemat- 
ical details of using the Empirical Probability Generating 
Function to provide estimators are found in section IV, along 
with a pictorial display of asymptotic efficiency for these 
estimators. An algorithm for computing the estimates of this 
method is presented in section V. The method of using the 
Empirical Probability Generating Function 1s compared to the 
Method of Moments in Section VI. The algorithm described in 
section V is illustrated by a numerical example along with 
computer studies in section VII. Tables for use with the 
above algorithm are found in section VIII and a summary of 


the results appears in Section IX. 





Il. NEYMAN TYPE A DISTRIBUTION 


Neyman's Type A distribution can be used to model decay 
and reproduction processes which fit the assumptions of com- 
pound Poisson processes. A derivation of this distribution 
provides an understanding of its structure, usefulness, and 
properties. For these reasons it is presented in the fol- 
lowing form. 

Consider a decay process in which radioactive nuclei 
decay to produce daughter nuclei according to a Poisson proc- 
ess with rate A’. Also suppose that the number of daughter 
nuclei produced after each decay are independent random 


variables having a common Poisson distribution with parameter 


m5. Let N(t) = the number of decays in the time interval 
Gone le: 
Y; = the number of daughters produced from the ith 
decay. 
Xx = the total number of first generation daughters 


produced in the time interval (0,t]. 


Then 
N (t) 
a y Y. 
i=l ~* 


Because N(t) is a Poisson process with rate i, 


=-At (At)” 
P IN(t) =n] = e n! I@ ee fg eke ea Oral ae ee rat 


ae 








Conditioning on N(t) is used now to find the probability mass 


function of the random variable X. For fixed N(t), X isa 


Poisson random variable with parameter nm. when N(t) =n. 
Thus 
x 
“nm, (nm. ) 
emo = x | N(t) =n) = e4 x! nor x = 1,25 35.508 


The event that no daughters are produced in the first gener- 
ation may be described by two cases, namely 
1) X= 0 if N(t) = 0 or, 
Ze Xk = 0 if N(t) # O. 
0 


From case 1) P(X = | N(t) = 0) = 1 is obvious. 


-nm., 
Poemecrrom case 2) P(x = 0 | N(t) =n#0) = e : 


Now, using the total probability theorem, 


Pee> P(X = 0) = & P(X = 0 Gi ce) a) 
n=0 
Or 
=it 0 “nm,  -At Oey 
P = e ~ ; e e n! 
= n=1 
= Tt 
-At(l -e a 
P = e 
Oo 


[ier also obvious that for X = x > 0, P(X = x | N(t) = 0) = 
And by using the same conditioning argument as above it is 


evident that, 





-nm. (nm Neos n 
00 2 2 “At (At) 


IDS 1S tes = oy e x! e De! i 
Xx 
hn 


These two formulas may be combined into one expression which 


is given in Shenton [3] as 


eas! Xx x x taped a. 
P = e m [0  eeele +ege2  +°8° 3° 8+ ct ee ] C2 ale) 
= z l mil 31 
x! ° ° 
where 
-m., Ih Shee eG 
my = At g = mje sleyel 0) = 

Oif x # 0 


This probability distribution was originally derived by 
Neyman and is called the Neyman Type A distribution. Other 
forms of equation (2.1) may be found in Johnson and Kotz [4]. 

Some of the properties of this distribution are needed in 
the further development of this paper and are 


1. Probability generating function, 


m., (t-1) 
Ke ee eee ] 
See. — (2.2) 
2. Moment generating function, 
m. (e° = 1) 
nx -m, [l-e ] 
E(e ) = e : (2353) 
By using the well-known formula 
ux 
k 
— Eee ao. = E(X*) 
du 


the following are obtained 





mym. (1 + m. + mm.) 


mm. (1 + 3m 


sy, 


m 


2 


Il 


2 
+ my mt. 


m, (1 +m 


+ 3m,m 


2 
=) 


>) 


10 


Jk 


2 


> 3m,m 


2 


2 


+ 


(2.4) 


(22) 


(276) 


(2.7 ) 





Iti. EFFICIENCY AND MAXIMUM LIKELIHOOD EQUATIONS 


A method for estimating the parameters of a distribution 
should yield estimators with high asymptotic efficiency and 
it is well known that Maximum Likelihood estimators have 
asymptotic efficiency equal to 1. However, the Maximum Like- 
lihood equations may be difficult to solve and such is the 
case for the Neyman Type A distribution. Replacing some of 
these equations in the Maximum Likelihood system may produce 
a system of estimating equations that is more readily solved 
and still produces estimators with high asymptotic efficiency. 
Estimating equations considered here are those that equate a 
Statistic with its expected value. This paper presents one 
such system for the Neyman Type A distribution. 

Tt can be shown (see, for example, Read [5]) that the 
asymptotic efficiency of the estimators of such a substitute 


system is 


-1 
err = |™ | 
| A 
where |A| = information determinant of the probability dis- 
tribution considered and 
“w “~~ Ji 
M = lim [ nE(@ - 8)(98 - 8) J] 
n-o 
with 9 = vector of distribution parameters 
and 98 = vector of estimators of 98. 


aE 





6 expression for EFF can be simplified if 






svector of estimating equations 


A = the matrix e (32) and C = lim [ nE (hh‘) ] 


n-7o 
Le 


Z 
[Al 
ee er 


JA] [c| 
esi mun Likelihood equations for the Neyman Type A 


ation are derived in [3] and the following forms are 


[4] 
x = Bes 
n Pp 
mee)|6=C (xX, «+ (1) Le Sia. i322) 
, a Pp 
| Soe 
a 


_= the sample mean 

= the sample size. 

Jshould be noted that equation 2) is an extremely com- 
me function of ms and m,- To compute the efficiency 
Other set of equations as substitutes for the equa- 
in (3.2), it is necessary to compute the information 
inant of the Neyman Type A distribution. The informa- 


Meeeminant has been derived in [3] and there it is 


that 


v2 





IV. ESTIMATING EQUATIONS 


The Maximum Likelihood equations are quite difficult to 
solve for the Neyman Type A distribution as seen from the 
equations in (3.2). Namely, the second equation is the one 
that produces the complication. Replacing the second equation 
in (3.2) with a different one will allow for simpler calcu- 
lations if this new equation is chosen judiciously. However, 
the idea of maintaining high efficiency should also be con- 
Sidered when making this choice of a substitute. The use of 
the Empirical Probability Generating Function (hereafter 
referred to as EPGF) seems to be appropriate in attaining 
this goal. This was suggested by the successful application 
of this choice in the case of the Negative Binomial distri- 
bution as demonstrated by Caglayan [7]. Since the EPGF 
contains the independent variable t, the solution of the new 
system of estimating equations can be "tuned" so as to achieve 
the estimators of maximum efficiency capable from the new 
system. Considering the form of the probability generating 
function for the Neyman Type A distribution, it is easily 
amendable to programming on current programmable calculators. 
With this motivation in mind, the first Maximum Likelihood 
equation in (3.2) and the EPGF equation which follows are 
considered as a possible alternative to estimating the pa- 


rameters for the Neyman Type A distribution. 
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The EPGF is defined as 


as 
1 


M5 


EPGF = 1 
nei 


where n = the sample size and x, = the ith sample value. The 
expected value of the EPGF is the probability generating func- 
tion of the distribution of the random variable X and is given 
by equation (2.2). The system of estimating equations con- 
Sidered above is then 


ax = m,m., (4.1) 


m. (t-1) 
x. “Mm, [l-e ] 


2) 2 =e (Ae) 


up 
te 


i 
n 


1=1 


The notation m, * and m.* is used to denote the eStimators 


of my and m. EOuUnGe by SOlyiIngseqmatvense(4.1) and (4.2)- |) fe 
solve these equations the variable t must have a value in the 
interval [0,1) to insure that the sign of the left-hand side 
of equation (4.2) is compatible with the right-hand side. 

The value chosen for t (which will be denoted as t*) is found 
by maximizing the asymptotic efficiency of the estimators m,* 


i 


and mo *. 


To compute the asymptotic efficiency of m, * and m. * 


equation (3.1) is used. To reduce the amount of typing in 
the following expressions let 
m. (t-1) 


-m, { lLl-e ] 
G(t) = e +t 
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The vector of equations 


Ay = 


The matrix A is 


(h) is given by 
x - mm. 
n Xs 

a8 ac - G(t) 
nee — 

dh, dh, 

dm dm. 

dh. dh. 

dm “| 


meemewhich it follows that its elements are 


2) ae 


oy = 


Ze 


Ie 


The matrix C is 


enennt 
n> 


and it follows that its 
Cia = lim nVAR(x) 
n+2 
C15 —eelim mCcOV (x, 
n+ 
the last step resulting 
ee 3. 


-m, CaS) 
-m, (Ad) 
l-e G(t) (aes) 
m. (t-1) 
m, (1-t)e @its) (4.6) 
2 
h, hh, 
oP UL eee cee 
Zo 2 
elements are 
= m,m., (1 + m. ) ao) 
n xX. x 
il Teo Po COV Xa 2 ont ) 
mae . L 
n i1=l 
rap azo) x and a being independent when 
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Semntinuing, 


Xx. Xx: 


Cc = E(x,t = ae B(x, )E(t +) 


1 
and using the fact that for well-behaved distributions (in- 


cluding Neyman's Type A) 


Xs x7 
E(x.t ") = tE [ a(t ~) ] = tG'(t) 
dt 
then 
m. (t-1) 
Ci = m,m.G(t) { te = le ae (4.8) 
By symmetry Ci2 = Coi- 
n Xs Xe 
C55 = lim nVAR(1l1 % t 7) = VAR(t ~) (4.9) 
n-© n i=l 
thus 
* 2 Z 
C55 = Git ) - [G(t)]~. (42) 


Using equations (4.3) - (4.6) 


m., (t-1) 
| = cutee + (mM. (cal) aalye ] (4.11) 
and using equations (4.7) ~- (4.10) with 
B = 1+ m. + m,m. (te = 1k) 
then 
[Cc] = mm, [ (1l#m,)G(t*) - G*(t) BJ. (4.12) 


Now substituting (4.11) and (4.12) into (3.1) results in 


m,G (eek + (mt - mM, - 1) e ] 
SPP = OOO Cees 


|Alm, { (1+m,)G(t*) - G*(t) Bi] 


ILy 





t* is the value of t used to solve equations (4.1) and (4.2). 
It is defined as 


EFF (t*) = max Siblae(Ge)) © 
Ose<k 


In order to use this definition the values of my and mM. 
must be given. t* and EFF(t*) have been computed for known 
mM and Mm. and are presented in Table 1. An algorithm for 
using these tabled values to find the estimates m, * and m. * 


is presented in section V.B. 


To describe the relationship between t*, mM, and M5, 


three-dimensional plot was made (Fig. 1) for mM and Mm. in the 


a 


mae. 01(.01).09, .1(.1).9, 1(1)9. Fig. 2 is a corresponding 
Prertrot EFF(t*) for the same set of values for my and M. - 
Figures 1 and 2 are representations of the values found in 
Table l. 

For given values of my and M5, PE tomenes value Of £ which 
Maximizes the efficiency function EFF(t). The values of t* 
are not necessarily continuous as my and m5 are varied. The 
efficiency function determined by equation (4.13) is bimodal 
for certain values of mM and m. and these modes occur near 
the endpoints of the interval [0,1). The global maximum of 
EFF(t) on this interval may occur at either mode for different 
values of m4 and Mo , which is the cause for the discontinui- 
ties of t* as seen in Fig. 1. The jump in the surface occurs 
when my 15 about 3.0 for all values of m, > 2.0. Further, it 
is interesting to note that t* is approximately zero for 


my < 1.0 and all values of Mo, except for a Singular point at 
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m, = ne) Ol ard m, = -O1. EFF(t*) appears to be a continuous 

mametlon of m4 and m5 as shown in Fig. 2 and is very close 

fom. 0 Over much of the range for my and mM. described above. 
The method used for determining t* is discussed in 


section V.A. 


lee 





O at the base 
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=) 
wm © 


parameters m5 and mM. 


im cc vs. 


akegr 
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V. CALCULATING TECHNIQUES 


Estimation of the parameters m4 and m- of the Neyman Type 
A distribution using the method outlined in section IV in- 
volves solving a nonlinear programming problem; specifically, 
finding a solution (m,*, m5*, t*) that maximizes equation 
(4.13) subject to the constraint equations (4.1) and (4.2). 
Since this requires a large amount of computing power, an 
approximating sequential method is proposed that can be per- 
formed on a calculator. In order to do this the values of t* 


must be known for any m, and m. pair. These t* values are 


il 
presented in Table 1 along with the associated maximum EFF 


for a wide range of values of m, and M. - Bivariate linear 


iL 
interpolation as discussed in Abramowitz and Stegun [8] can 
be used to approximate t* between the values of mM, and m. 


given in Table l. 


A. METHOD OF COMPUTING t* FOR TABLE Il 

Table 1 was prepared by fixing m4 and Mm. at constant 
values in the range .01 to 9.0 and then conducting a single 
variable search on the function EFF (as a function of t only). 
A Golden section search procedure was employed that terminated 
at the value of t* equal to the midpoint of the searched in- 
terval on which the value of EFF(t) was increasing but did not 


change by more than 107 across the interval. Because of the 


bimodality of EFF(t) on some regions, a search was initiated 


ZZ 





from the right end of the interval [0,1) and a separate search 
started from the left end. The global maximum was then found 


by comparing the resulting local maximums. 


B. METHOD OF SOLVING FOR m, * AND m. * 

Given a sample of n observations on the random variable X 
assumed to have an underlying Neyman Type A distribution of 
which x, represents the ith sample observation, the algorithm 
for solving equations (4.1) and (4.2) to find the estimators 
m, * 


lowing steps: 


and m. * of the parameters m, and m5 consists of the fol- 


STEP 1 - Compute the Method of Moments estimators from 


the following equations 


_ n 
x = a5 > x, 
n i=l 
n 
See ae Go Hc isa 
n-1l 1=1 
m, = x/M. 
m S29 
m, = =a al 
X 


STEP 2 - Find t* in Table 1 corresponding to m4 and M. - 
ieee If t* = 0. and the sample contains no zero values, set 
m,* = my and m,* = m. then STOP. 

STEP 3 - Solve the system of equations (4.1) and (4.2) 
ror m, * and m. * with t = t* by using a Newton-Raphson pro- 
cedure. Solving equation (4.1) for my and substituting this 


into equation (4.2) yields 


Zs 





where 
: m. (t-1) 
Gg = e : 
Applying the Newton-Raphson method on equation (5.2) produces 


as the (k+l)st iterate of m.* 


m, (t-1L) — = 
1l1-e K - my. In (t*) 
1+e (mt -™m - i) 


xX 
m = m [1+ 
2 (k+1) K 


where for ease of typing nm. = m and m, = m., (the Method 


2 
(k) Cis) 

Of Moments estimator). When termination of the above iteration 
scheme occurs by the appropriate criteria as set by the user, 


set m. * = m, and compute m, * Prom eGuataem (5) 


(n) 

STEP 4 - If |m.* - Mm. | Qe STOP wheteue mic cere by atac 
user. Otherwise return to STEP 2 after replacing m5 with m,*. 
The convergence of this algorithm depends on the sample 
Size n and how well the Method of Moments eStimators approxi- 
mate the true values of mM, and M. Experience has shown that 

cycling between estimates can occur when the m,* and m.* 
values are close to the points where t* is discontinuous as 
determined from Fig. 1. This behavior was observed when the 
above algorithm was applied to data of sample sizes as large 


as 50. For most cases the estimates obtained by applying the 
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Megeorithm converge after two or three iterations. It should 
be noted from STEP 2 that if t* = 0 and the sample contains 
no zero values, then the Method of Moments estimators are the 


default values obtained when using the algorithm. 


feet HOD OF DETERMINING Q 
The value of Q as defined in equation (3.4) was computed 

eacecording to equation (3.5) for the purpose of determining the 
asymptotic efficiency of the EPGF estimators given in Table l. 
The xth power moment of a Poisson distributed random variable 
(denoted by Si in equation (3.5)) was computed by finding the 
maximum term of the series and using this as a scaling factor 
to keep the partial sums within the limit of an IBM 360/67 
computer, which was used for all computations. Actual compu- 
tations were performed with 16 digit precision using logarithms. 
Truncation of the series occurred when the first scaled term 


~lel.l. The 


after the maximum term was found to be less than e 
value of Q was calculated by computing the successive partial 
sums of ae. until the relative error between the last two was 


7 


less than 10 ‘ (the last term was not added to the partial sum 


computed previous to it), 1.e., when 


Im+1 _7 
< 10 ; 


Peetiseussion of this stopping rule for the particular appli- 
Cation described above may be found in Katti and Gurland [9]. 


As seen in Table 1 there are a few instances where the value 


ZS 





of EFF (t*) was computed as a value greater than 1.0,which is 
impossible. It is felt by this author that in these occur- 
rences the value of Q was not computed accurately enough. 

The reason for this is that when m, Or m, was large over 300 
terms were required in the partial sum of Q in order to reach 
the stopping rule criterion given above. However, since the 
size and frequency of these inconsistencies were small, it is 
felt that the other results presented in Table 1 are not com- 
promised. This is justified by noting the continuous behavior 
of EFF(t) as given by equation (4.13) and observing the few 
Values of EFF(t*) > 1.0 in Fig. 2. It is not known why these 
occurrences did not happen at the extreme values of m, = O50 
and m5. = .0O1. As shown in Fig. 2, they appear to be isolated 


Semmes OL irregularity in the computation of EFF(t). 
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VI. COMPARISON WITH ESTIMATION BY THE METHOD OF MOMENTS 


To determine the utility of the EPGF method of estimation 
for the Neyman Type A distribution, a comparison with Method 
of Moments estimation was made. This comparison was performed 
by computing the ratio of the asymptotic efficiencies of the 
two schemes. The estimators of my and m, by the Method of 
Moments are found by solving the equations of (5.2) and is 
actually done in STEP 1 of the sequential algorithm presented 
in section V.B. The EPGF method is an attempt to improve on 
the Method of Moments estimation of M4 and M5 - The amount of 
improvement is determined by the increase in asymptotic ef- 
ficiency gained by applying the EPGF method relative to that 
of the Method of Moments. In order to find this increase, 
the asymptotic efficiency of the Method of Moments must be 
found for the Neyman Type A distribution. Using equation 
(3.1), this was done and is shown below. 

The determinants of A and C for the Method of Moments are 
determined in a manner similar to section IV. By using the 


equations in (5.2) the following results can be obtained: 


Shy areal 

a en 

oe ee ma” 
A, 5 = ~m, ~ 2m,m. 


Zi 





thus 


_ 2 
|A| = mm." . Meek) 
Also 
Chi = VAR (X, ) = mm. (1 7 m.) 
Oe E(X.°) - 36(x.°) E(x.) + 2 [ B(x.) 3° 
i an a 5 1 
c = mi,m,(l + 3m, + m 2) 
2 i 2 Z 2 
and 
te %12° 
From formulas derived in [5], 
pee (x. —- X)° - { vaR(x.) }¢ 
22 at af 
C =" mm, { 2m.m,.(1 + m \ 2 + (1 + 7m. + 6m 2 + m 3) ] 
22 gee 2 2 Z 2 2 
so that 
fy Ze 2 Jee 2 3 
fe) = Mm, | 2m,m4(1 +m) + 2m, + 2mo° +m” J. (6.2) 


Substituting equations (6.1) and (6.2) into equation (3.1) 
yields 
Z 
a 
BEF Ss ST on) 


> 
JA[ { 2m,m,(1 + m,)~ + 2m, + 2m," + m 


where EFF am represents asymptotic efficiency of the Method of 
Moments estimators m and M. - The ratio EFF/EFF formed from 


equations (4.13) and (6.3) represents the relative merit of 


the EPGF to Method of Moments estimation. 
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Values of t* and EFF/EFF are tabulated in Table 2 for 
the same range as covered by Table 1. This comparison of 
the two methods is most readily observed in Fig. 3, which is 
a three-dimensional representation of Table 2. The EPGF 
method is better than the Method of Moments in producing 
estimators of high asymptotic efficiency for the range of 
parameter values covered in Table 2. This is clearly demon- 
strated by Fig. 3. The EPGF method reduces to the method of 
estimation by zero frequency when t* = 0. In this case the 


contours of Fig. 3 can be found in [4] and [9]. 
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bee ( t * ) 
ayaa 
Tum 


Fig. 3 Efficiency of the EPGF method relative 
to the Method of Moments vs. the 
parameters m4 and M5 
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VII. NUMERICAL EXAMPLES 


EXAMPLE 1 

To illustrate the EPGF method and in particular the al- 
gorithm of section V.B, the data of insurance claims caused 
by accident or sickness cited by Shenton [3] are used and 
reproduced below. The data are assumed to be observations 
from a Neyman Type A distribution and three methods of esti- 
mating the parameters m5 and m5 are considered; 

1) the Method of Moments, 

2) the EPGF method, and 

3) the Method of Maximum Likelihood. 


The data are: 


Value Frequency 
0 Vey 
1 lees Ss 
2 200 
5 164 
4 O77 
5: 68 
6 49 
Wi 39 
8 2 
9 V2 
10 Gab 
itak 2 
We 5 
3 Z 
14 3 
is ul 
sample size = 1056 
sample mean = 2200 59 
sample variance = 6.4106. 
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Paoneede EPGE algorithm of section V:.B is applied to this 
data, the following results are obtained: 
a) from STEP 1 the Method of Moments estimators are 
computed to be 
mM = 2.184 and m. = 1.285 


Bb) from STEP 2 the t* value corresponding to my and Mm. 


is found by bivariate linear interpolation in Table 1 as 


follows; 

for m, = 2.0 and m5 = i2o0) Linear ancerpolatlonedives 
t* = .0748, 

Or m, = a.0 and m, = 1.285 linear interpolation gives 
mo 2.35/74, and 


1a o hg m, = 2.184 and m. =1.285 linear interpolation gives 
foe .1305 

“meeectrom STEP 3 the solution of (5.2) yields 

m,* = 2.638 and m. * = 1.064 

where termination of the Newton-Raphson method was set for an 
Seeerice error of .0001 

da) from STEP 4, where ¢€ = .01, a second application of 
the algorithm is required and the above m,* and mM. * 1s used 
fememe= new starting point in STEP 1. This second iteration 
then yields 

Eo 2736, m, * =) 2609. alc m. * = 1.076. 

Since the termination criteria of STEP 4 is not met a third 
iteration is performed resulting in the following 


cao. 2093, ~ = 22,611, and m. * = 1.074 


ma 


with the termination criteria being satisfied. 
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The computation of Maximum Likelihood estimators for the 
Neyman Type A distribution is very complicated and only ap- 
proximate Maximum Likelihood estimates for the above data 
are given in [3]. To get an idea of the computational effort 
involved, these estimators are only the first iteration re- 
Sults and required several hours on a desk calculator in 
}949. 

Maximum Likelihood estimators for the above data were 
found using a TI-59 programmable calculator and the method 
presented in [6]. The range of observed values was small 
enough so that data storage was not exceeded on the calculator. 
The following results were obtained in approximately thirty 


minutes of computing time, 


Tteration ak 2 
ne 2a 5 1.134 
Z 22515 jiu als 
3 2 ke deel 
4 or ako i lapes as L 


In comparison, the computing time used to perform the 
EPGF method was approximately five minutes and the data 
storage requirement consisted of only keeping eight numbers, 


which were not influenced by the data values. 


EXAMPLE 2 

Computer simulation was used to generate eighty-one data 
sets of random numbers from the Neyman Type A distribution. 
Computer program 4 in Appendix A was written to simulate the 


model presented in section II and the data sets produced are 


Sh. 





given in Appendix B. This program uses a random number gen- 
erator that provides random deviates from a Poisson distri- 
bution and is a local subroutine. For each data set produced, 
the sample size was fifty and the range of values used for mM 
and m, was TE) ae a 

The EPGF estimates for each data set were found by apply- 
ing the algorithm presented in section V.B. The t* value 
required in STEP 2 was computed by actually maximizing the 
efficiency equation (4.13) with the estimated values of mM 
and M, « The results of these computations are found in Table 
3. The termination criteria of STEP 4 was set at « = .0001, 
so that more iterations were required than would normally be 
expected. 

In only two of the eighty-one cases considered, no im- 
provement was achieved over the Method of Moments estimators. 
These two occurrences are noted in Table 3. For the case 
m, = 4.0 and m, = 7.0, the special termination of STEP 2 of 
the algorithm was invoked; i.e., t* = 0 and the sample con- 


tained no zero values. For the case m, = 5.0 and m5. = 650), 


cycling between the following estimate values occurred, 








tx m,* mo * 
.0000 Bei 9.6798 
nomi o 2.9465 10.5754. 


The successive estimate pairs m,* and m. * produce values 
of t* which are on either side of the discontinuous ridge 
seen in Fig. 1 of section IV. As previously discussed, the 


efficiency function is bimodal in this region and causes the 
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above behavior. Since the asymptotic efficiency is high in 
this region, a larger sample size producing better initial 
estimates should resolve the inconsistency which results 


when cycling occurs. 
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VIII. TABLES 


Tables 1, 2 and 3 mentioned in the previous sections are 


contained on the following pages. Table 1 presents the op- 


timal t and efficiency values found for the range of my and 
m. from .01 - .09 (increments of .01), .1- .9 (increments 
O@ie.l), and 1.0 - 9.0 (increments of 1.0). Table 2 presents 


the optimal t and efficiency (relative to the efficiency of 
the Method of Moments) values found for the same range as 
used in Table 1. Table 3 presents the results of applying 


the algorithm of section V.B to the data sets of Appendix B. 
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OPTIMAL T ANDO EFFICIENCY VALUES 


TABLE 1: 


M2 VALUES 
0.0500 


Ml VALUES 


0.0700 0.0800 0.0900 


0. 0600 


O.-C2CC 0.0300 0.04C0 


0.0100 


0.01 


0.02 


0.04 


0.05 


0.06 


oy) 


0.07 


om 
oO” 
Oo 
oOo 
e © 


oOo 


oo 
Oo 
om 
OO 


oO 


© O° 


OOo 
oOo 


eo @ 


0 
0 


0. 0000 
0.9999 
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0.9999 


O01 
0000 


One 


mO 
i 
oq 
ore) 


NO 
Oo 
Tan) 
20 


C 
1 


0.0003 
1.0000 


0.08 


Ct) ad 
La Cr) 
DP mt 
SO 

e@.hU6°8 


Oa 


0.09 


ENTRY IN THE TABLE: 


FOR EACH 


EACH PAIR 


TOP VALUE OF 
BOTTCM VALUE OF EACH PAIR 


J 


EFF 





OPTIMAL T AND EFFICIENCY VALUES 


FABLE 1: 


M2 VALUES 
0.5000 


Ml VALUES 


0.7000 0.8000 0.9000 


0.6000 


0-2-2000 0.3000 0.4000 


0.1000 


0.01 


0.02 


0.0000 
0.9998 
0.0000 
0.9997 


oOo 
oOo 
Oo 
Oo 


0 
0 


0.0000 
0.9999 


0.03 


ead CD 
oo 
Oo 
Ooo 

ee 


or 


0.0% 


0.05 


0.06 
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0.0000 
0.9999 


0.07 


or 
OF 
oOo 
Oo 
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Oo 
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aloo) 
om 
on 

e @ 


Ooo 


ovr 
ao T so) 


oo 


OO 


Ot 
On 
Oo 
oOo 
e 8 


O 
0 


0.0000 
0.9990 


0.-00CO 
0.9993 


0.02 


om 


0.09 


FOR EACH ENTRY IN THE TABLE 


TOP VALUE OF EACH PATR 


Tx = 


BOTTCM VALUE OF EACH PAIR 


— 
—= 
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OPTIMAL T AND EFFICIENCY VALUES 


TABLE 1: 


M2 VALUES 
5- 0000 


ML VALUES 


#.0000 8.0000 9.0000 


6.0000 


2-0C00 3.0000 4.0000 


O 
0 


1.0000 
0.0000 
0.9999 


0.01 


0.02 


0.03 


Ce 
On 
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0.04 


0.05 


0.0000 
0.9982 


0.06 
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0.0/ 


0.08 


0.09 


FOR EACH ENTRY IN THE TABLE 


TOP VALUE QF EACH PAIR 


Tx 


BOTTCP VALUE OF EACH PAIR 


EFF 





OPTIMAL T AND EFFICIENCY VALUES 


TABLE 1: 


M2 VALUES 
0.0500 


Ml VALUES 


0.0700 0.0800 0.0900 


0.0600 


0.0400 


0.0300 


Q-O0L00 0.0200 


0.10 


0.20 


0.30 


0.40 


0.50 


0.60 


40 
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0.90 
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VALUES 


OPTIMAL T AND EFFICIENCY 


TABLE 1: 


M2 VALUES 
0.5000 


Ml VALUES 


0.7000 0.8000 0.9000 


0.6000 


Ce-2000 0.3000 0.4000 


0- 1000 


0.10 


0-20 


0.30 


0.40 


0.50 


0.60 
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OPTIMAL T AND EFFICIENCY VALUES 


TABLE 1: 


M2 VALUES 
5.0000 


Ml VALUES 


7.0000 8.0000 9.0000 


6.0000 


2-C000 3.0000 4.0000 


1 -Q000 


0.10 
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o~ 
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0 


0 0000 
0.9862 
0.0000 
0.9866 


0.50 


ON 
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0 
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0.90 


FOR EACH ENTRY IN THE TABLE: 


TOF VALUE OF EACH PAIR 


Ts 
EFF 


OF EACH PAIR 


BOTTOM VALUE 





OPTIMAL T AND EFFICIENCY VALUES 


TABLE 1: 
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FOR EACH ENTRY IN THE TABLE: 
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BOTTCM VALUE OF EACH PAIR 


— 
—_ 


EFF 





OPTIMAL T AND EFFICIENCY VALUES 


TABLE 1 


M2 VALUES 


M1 VALUES 


9.0000 


720000 8.0000 


2-2-0000 3.0000 4.0000 5.0000 6.0000 


1 -0000 


1.00 
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3-00 
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FOR EACH ENTRY IN THE TABLE 
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EFFICIENCY VALUES RELATIVE TO THE METHOD OF MOMENTS 


OPTIMAL T AND 


TABLE 2: 


M2 VALUES 
0.0500 


M1 VALUES 


0.0900 


0.0700 0.0800 


0. 0600 


C€.0200 0.0300 0.04C0 


0.0100 


0.01 


Om 


0.02 


0.03 


oF 


0.04 


0.06 
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0.07 
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0.09 


FOR EACH ENTRY IN THE TABLE: 


TOP VALUE OF EACH PAIR 


Tx 
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BOTTCM VALUE OF EACH PAIR 





THE METHOD OF MOMENTS 


GPTIMAL T AND EFFICIENCY VALUES RELATIVE TO 


TABLE 2: 


M2 VALUES 
0.5000 


Ml VALUES 


0.7000 0.8000 0.9000 


0.6000 


0.2000 0.3000 0.4000 


0.1000 


0.01 
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FOR EACH ENTRY IN THE TABLE: 
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OPTIMAL T AND EfFICTENCY VALUES RELATIVE TO THE METHOD OCF MOMENTS 


TABLE 2 


M2 VALUES 


5.0000 
0.0000 
1.5384 


M1 VALUES 


7.0000 8.0000 9.0000 


6.0000 
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OPTIMAL T AND EFFICIENCY VALUES RELATIVE TO THE METHOD OF MOMENTS 


TABLE 2 


M2 VALUES 


M1 VALUES 
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IX. CONCLUSIONS 


The EPGF method provides estimators with asymptotic ef- 
ficiency close to that of the Maximum Likelihood estimators 
_over most of the parameter values considered. The smallest 
value found was 54% (at m, = 3.0 and m5, = 9.0) and in most 
instances the asymptotic efficiency is greater than 97%. 

The EPGF method outperforms the Method of Moments in pro- 
ducing estimators of consistently higher asymptotic efficiency 
for all parameter values from .01 to 9.0. 

The efficiency of the EPGF method relative to the Method 

of Moments is less than 110% when m < -3 for most values of 


My - As may be seen from Table 3, the relative efficiency 


drops below this value again for large my and m5 beginning 


about m, = 9.0 and m, = 3.0, and extending to m, = 6.0 and 
Mm. = 9.0. 
In almost all cases, t* = 0. when my is less than 1.0. 


When this occurs the EPGF method reduces to eStimation by the 
Method of Zero Frequency. The EPGF method becomes an exten- 
Sion of the Zero Frequency method when the efficiency of the 
estimators of the latter decreases. This happens when t* > 
0. which is always the case for m, > 3.08 

Although the EPGF method is always better than the Method 
of Zero Frequency in terms of efficiency, it involves a 


greater computational effort when t* # 0. As shown in [4], 
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the Method of Moments estimators have higher asymptotic ef- 
ficiency than those of the Method of Zero Frequency whenever 
m, > B> Lor all values of M5 - The efficiency of the Method 
of Moments estimators is over 20% higher than the Zero Fre- 
quency estimators when m, > 4.5. Combined with the above 
discussion of the EPGF method, whenever t* > 0. the EPGF 
method will yield estimators of substantially higher asymp- 
totic efficiency than those found by the Zero Frequency method. 
The choice of which method to use always rests with the user. 
However, if it is assumed that m, > 22) aligns! wavehe Jee e)9 m4 and 
m5 are not large (greater than about 7.0) then the extra ef- 
fort required by the EPGF method is compensated for by a 
better than 10% increase in asymptotic efficiency of the 


estimators. When it is supposed that m, and m. are not to be 


iL 
found in the above region then estimation by the Method of 
Moments is recommended, as little increase in efficiency can 
be expected by using the EPGF method; i.e., going beyond 
SHEP 1 of the algorithm presented in section V.B. 

As noted in section VII, the EPGF method as applied by 
the algorithm of section V.B can fail to converge. fThis is 
characterized by cycling between two different pairs of 
estimators when the sample size is small. When this occurs 
the remedy is to use a larger sample. The variance of a 
population having the Neyman Type A as an underlying distri- 


bution grows large as m4 and m. increase, especially when 


both are greater than 1.0. Assuming that the sample size is 


oul 








large enough for the Central Limit theorem to provide a good 
approximation by a Normal distribution, the sample size re- 
quired to provide a reasonable confidence interval containing 
the true mean is quite large. Since x = mm, is one of the 
estimating equations used in the EPGF method, it may be dif- 
ficult to provide good initial estimates for the iterative 
procedure when n is small. It has been found by experience, 
as the sample size iS increased, the EPGF estimates tend to 
cluster more tightly than the Method of Moments estimators 
about the true parameter values. 

The EPGF method is not as easy to use as the Method of 
Moments, but with the aid of Table 2 and the use of the 
algorithm presented it can be applied readily with current 
programmable calculators. This is not always the case for 
estimation by Maximum Likelihood. 

Although the EPGF method appears to be time consuming 
and tedious, it is far easier to implement than the Maximum 
Likelihood method. By using Table 1, reasonably good esti- 
mators can be obtained with less effort than required by the 


latter method. 
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APPENDIX A 
COMPUTER PROGRAMS 


Four computer programS were used in the compilation of 
this paper and are given on the following pages. All were 
programmed in the standard Fortran IV language and in double 
precision. Program 1 was used to compute the values of the 
information determinant |A| of section III for m, and m, in 
the range of Table 1. Program 2 was run to compute the 
values of t* and EFF(t*) and produce Table 1. Table 2 was 
produced by slightly modifying program 2 to compute and table 
the efficiency ratio EFF/EFF on of section VI. Program 3 was 
used to perform the algorithm of section V.B and produce the 
estimates for the numerical examples of section VII, which 
are found in Table 3. Program 4 was used in Simulating ran- 
dom numbers from the Neyman Type A distribution. The sub- 
routine LPOIS1 found in program 4 1S a local random number 
generator for Poisson distributed random deviates. 

Programs 1, 2, and 3 are complete and may be used di- 


rectly. Program 4 must be modified to be compatible with 


the user's local random number generation subroutines. 





COMPUTER PROGRAM 1 


G COMPUTATION OF INFORMATION DETERMINANT GIVEN Ml AND M2 
Pie bbe tr REAL“ (A-H,;O-Z) 

REAL*8 M1,M2 

CeMMON TMI M2, BETA, DET 

DIMENSION PARM(27) ,DETQ (27,27) 

Pai PARM/ beD—-2,2.D=2,3.D-2,4.D=2,5.D=2,6.D-2,7.D=2, 

Pee ee, Dt 2 P=. D1, 4.D=1,5.D=-1,6.D=-1,7.D-1, 
ecm oD 1 7) .DO,2. DO, orb0,4.D0,5-00,6.D0,7 .D0 8.00, 
Peo O / 

CREATE A TABLE OF DETERMINANT VALUES FOR VARYING 
PARAMETERS OF 

Ml AND M2 


QAAN 


M2 = PARM(J) 
BETA = ML*DEXP (-M2) 
VALUES HAVE BEEN SET FOR Ml AND M2 


AANA 


COMPUTE INFORMATION DETERMINANT 
CALL DETOS 
DETOcl,J) = DET 
99 CONTINUE 
feniteNe. 1) GO TO 150 

Neer (6,101) (PARM(L) ,L=J1,32) 
0 1 MORMAT (5X,3(3X,F12.10)) 
150 WEEE (6,151) ML, (DETO(L;L), L=Jig7s2) 
eS MBORMAT(F5.2,3(3X,F12.10) ) 
999 CONTINUE 
STOP 
END 


SUBROUTINE DETQS 
PMP LICL T REAL*8 (A~H ,O-Z) 
REAL*8 M1,M2 
COMMON M1,M2,BETA,DET 
C COMPUTE Q GIVEN Ml AND M2 
i =e. b= 7 
Po Pex (ae = ML) 
RX BETA 
QXx SRD IES 5k 
Se, oA 
ELN = DLOG(BETA) + BETA 
PO= Ogrkx) = 17500 
X = DFLOAT (IX) 
ELNEXT = SERIES (X+1.D0) 
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10 
20 


me 
Io 


20 


30 


40 


30 


RNEXT = ESXP(ELNEXT - ELN) 
ELN = BLNEXT 

QX = QX*M2* (RNEXT/X) * (RNEXT/RX) 
RX = RNEXT 

ieox/7O .bT. H) GO TO 20 

K lp, 

Ome Ox 

M2*M2*PO#Q 

M1L*M2*M2 *M2 

M1*M2*M2* (M1+M2+M1*M2) 
DET = (1.D0 + M2)*Q/A - B/A 
RETURN 

END 


Q 
Q 
A 
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Bune TLON SERIES (XxX) 
Pie iteit REAL*8 (A-H,0O-2) 
REAL*8 M1,M2 
COMMON Mi>M2,BETA,DET 
SERIES = 0.D0 
FIND INDEX OF THE MAX TERM: AIM 


AIM = 1.D0 
Gis= 1.D0 
e2y— 2.D0 


beelvcd K = 1,100 

FNUM = C2*AIM*DLOG(BETA/AIM) + C2*X - Cl 
FDEN = C2*(X + AIM) - Cl 

AINEXT = AIM*(Cl + FNUM/FDEN) 

[ieabNexe . Le. 1.D—5)° GO TO lil 

AIM = AINEXT 

CONTINUE 

IF(AIM .LE. 1.D0) AIM 
PeAIM ]GrT. 1.D0) AIM 
fr = [DINT (AIM) 
COMPUTE MAX TERM OF THE SERIES 

ALNIM = 0.DO 

Boe2o Tf = 1,IM 

Bae = DE LOAT (1) 

ALNIM = ALNIM + DLOG(AI) 

ALNIM = X*DLOG(AIM) + AIM*DLOG(BETA) - ALNIM 

COMPUTE EACH TERM (K) OF SERIES AND SUM TO THE SERIES 
bowlo ne— 1,300 

AK = DFLOAT (K) 

BENE = 0.D0 

PencOs = 1,K 

ie DELOAT (J) 

ALNK = ALNK + DLOG(AJ) 

ALNK = X*DLOG(AK) + AK*DLOG(BETA) - ALNK 

BK = ALNK - ALNIM 

IF (BK .GE. -161.1D0) SERIES = SERIES + DEXP (BK) 

Meni. Et a —1oi.1DOAND, K .GT. IM) GO TO 50 
CONTINUE 

COMPUTE THE NATURAL LOG OF THE XTH MOMENT. 

SERIES = ALNIM + DLOG(SERIES) 

RETURN 

END 


126 
DELOAL Clon ei 0 |) 


Wool 


65 





AMT Aa a1 OQ 


50 
HPO 
ei 
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COMPUTER PROGRAM 2 


COMPUTATION OF OPTIMAL T AND EFFICIENCY GIVEN Ml AND M2 
Meer REAL 8 (A—-H,O-—Z) 

REAL*8 M1,M2 

COMMON M1,M2,DET,TSTAR, ESTAR 

BOMENSTON PARM(27),TOPT (27,27) ,EOPT (27,27) ,DETQ(27 ,27) 

CREATE A TABLE OF T-OPTIMAL VALUES FOR VARYING 
PARAMETERS OF 

Ml AND M2 


READ IN Ml AND M2 AND CORRESPONDING INFORMATION 
DETERMINANT VALUES FROM A DISK 


i=) 0 

DO 99 K = 1,9 
Jl = J2 +1 
2 = ~J) + 2 
Deego © = 1,27 


Hee NE. 1) GO TO 50 

mew (5,101) (PARM(J) ,J=J1,52) 
eee oe VOL) (DETO(I,J) ,J=d1,02) 
BemMayT (5X,3(3X,F12.10)) 

CONTINUE 

SET Ml AND M2 VALUES FROM PARM 
COMPUTE OPTIMAL T-STAR = TOPT AND OPTIMAL EFFICIENCY 
E-STAR = EOPT 

Wer no9 Tf = 1527 

Ml = PARM(T) 

DO 199 J = 1,27 

M2 = PARM(J) 

mer = DETO(L,J) 

CALL OPTIML 


TOPT(I,J) = TSTAR 
EOPT (I,J) = ESTAR 
CONTINUE 


GeiNt OUT TABLE OF OPTIMARSY AND EEPECIENCY 
IU = Q 


DO 299 KI = 1,3 
IL = IU + 1 

IU = IL + 8 

JU = 0 

DO 298 KJ = 1,3 
JL = JU+i1 

JU = JL + 8 


Wee (6,20) 1) 

WRITE (6,202) 

WRITE (6,203) 

WRITE (6,204) (PARM(J) ,J=JL,JU) 

DOmezeo, L = 1L,10 

WRITE (6,205) PARM(I), (TOPT(I,J) ,J=JL,JU) 
WieteeG, 206)" (HOPN(L ad) 3-0 L,00) 
CONTINUE 
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WRITE (6,207) 
WRITE (6,208) 
WRITE (6,209) 

298 CONTINUE 

299 CONTINUE 

201 FORMAT('1',//////) 

202 FORMAT('0O',33X,'TABLE 1: OPTIMAL T AND EFFICIENCY 
* VALUES',//) 

203 FORMAT('0',9X,'Ml VALUES',35X,'M2 VALUES') 

204 FORMAT('0',19X,9(2X,F6.4)) 

205 FORMAT('O',11X,F4.2,4X,9(2X,F6.4)) 

206 FORMAT(' ',19X,9(2X,F6.4)) 

207 FORMAT('O',/,22X,'FOR EACH ENTRY IN THE TABLE:') 

208 FORMAT('0',21xX,'T* = TOP VALUE OF EACH PAIR') 

209 FORMAT ('0O'!,21X, 'EFF = BOTTOM VALUE OF EACH PAIR‘) 
STOP 
END 





SUBROUTINE OPTIML 
IMPLICIT REAL*8 (A-H,O-Z) 
REAL*8 M1,M2 
COMMON M1,M2,DET,TSTAR,ESTAR 
BEAGE= "0 
THETA = 9999.D-4 
Bee 0.DO 
C FIND A COARSE BRACKET AROUND THETA: (A,B). 
5 THETA = THETA - 1.D-2 
PATHE tA. LT. 0.D0) GO TO 15 
ET = EFF (THETA) 
fwibueehT, £2) GO TO 15 
H2 = ET 
eo. Tr 5 
5 fe) THETA 
fama rc O- DO) A = 0.D0 
B = THETA + 2.D-2 
IF(B .GT. 9999.D-4) B = 9999.D-4 
C 
PERFORM A GOLDEN SECTION SEARCH FOR TSTAR TO MAX 
C EFF (TSTAR) 
SeewiTH STARTING INTERVAL (A,B). 
R= (-1.D0 + DSQRT(5.D0))/2.D0 
20 Al (1.D0 - R)*(B - A) +A 
Bl R*(B - A) +A 
EA1L = EFF (AL) 
EB1L = EFF(BL) 
IF (DABS(EAl - EBl) .LT. 1.D-7) GO TO 30 


Powe ecr. Bl) B= ei 
TF(EAL .LT. EBl) A = Al 
GO TO 20 
30 Peewee Ae Bi) 2 bo 
ESTAR = EBL 
Cc CHECK FOR A MAX FROM THE LEFT END 


IF(FLAG .EQ. 1) GO TO 40 
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55 


36 


40 


he (TSTARY LT. 12D RETURN 


T = TSTAR 
E = ESTAR 
THETA = 0.DO 
E2 = Q.D0 


THETA = THETA + 1.D<2 

IF (THETA .GT. 99.D-2) RETURN 
ET = EFF (THETA) 

Magee, <n. EZ) GO TO 36 

b2e= ET 

SOe5O 35 

Pen elLubTA = 2.D=2 

tA aut. 0.D0) A = 0.DO 


B = THETA 

FLAG = lL 

SOerO 20 

Tete 2. LP. ESTAR) RETURN 
TSTAR = T 

ESTAR = E 

RETURN 

END 


FUNCTION EFF (T) 
BMP nL ellT REAL*~8 (A-H,0-Z) 
ReAGe Ss MI ,M2 
SeMMon MI,M2,DET,TSTAR, ESTAR 
COMPUTE PROBABILITY GENERATING FUNCTION G(T) 
pe DEXP (M2*T = M2) 
See OEAr (—-Mil*(1.D0 — F)® 
SOMPUTE G(T**2) 
hee DEXP (M2*T*T — M2) 
cee —e DE xP (—-Mi*(1.D0 = F2)) 
COMPUTE NUMERATOR OF EFF (T) 
ie |. bOet (M2*T — M2 — P2p0) sr 
Bie= ML *G*G*FN*EN 
COMPUTE DENOMINATOR OF EFF (T) 


Pt DO aM. ate Me Ae epee CLP = 1. DO) 
BDe= M2*((1.D0 + M2) *G2Z2 = G*G*FD) 
BDe= DETED 


COMPUTE EFFICIENCY EFF (T) 
EFF = FN/FD 

RETURN 

END 
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COMPUTER PROGRAM 3 


ce ESTIMATION OF Ml AND M2 


IMPLICIT REAL*8 (A-H,O-Z) 
DIMENSION 1X(450,10) ,M2(9) 
COMMON HATM1,HATM2 
Ree (59) (M2(5) ~=1, 9) 
9 FORMAT (3X,915) 
Cc READ IN MATRIX OF RANDOM NOS.: IX 
Reems, bo) ( (IX(1,3) ,g=1, 10), 1=1,450) 
10 FORMAT (I3,9IS) 
NNN 50 
DF LOAT (NNN) 
IFLAG = -l 
Pong) - = 1,9 
IFLAG = IFLAG* (-1) 
IK = NNN*¥(I-1) + 1 
IK + NNN - l 
PM1 = DFLOAT(IX(IK,1)) 
foeEP GAG .LT. 0) GO TO 15 
WRITE (6,100) 
100 FORMAT ('1',////,T37,'TABLE 3: NUMERICAL RESULTS') 
WeEte (6,101) 
101 Pee Ole / To MIN 127 No Tes Mi. T4353 5'M2" 151, 
fen Poe. ' Te! T65 Mie 7 oo 2+ reo. NO. ) 
WRITE (6,102) 
iL RORMAime ster sd a(MM)" TA2s (MMe oO aN’ )}T56,'FINAL', 
aersg; TER. ) 
5 CONTINUE 
BOowlo. — 2,10 
w= Fb — 1 
PM2 = DFLOAT (M2 (JK) ) 
XBAR = 0.DO 
MSO = 90, DO 
BO.20 K = IK,IKE 


re 

A 

3) 
1 


XBAR = XBAR + DFLOAT (IX (K,J) ) 
20 mo@1— xo) + DFLOAT (IXIK, J) *ixX(kK, 3) 
XBAR = XBAR/ANN 
XVAR = (XSQ - ANN*XBAR*XBAR)/(ANN - 1.D0) 


HATM2 = (XVAR - XBAR) /XBAR 
HATM1 = XBAR/HATM2 


AM1 = PML 
AM2 = PM2 
BM1 = HATM1 
BM2 = HATM2 
i PiND ESTIMATORS USiNG EMPIRICAL PGF 
PAST = HATM2 
JJJ = 0 
30 CALL TOPT (T) 


Peo 2 Ore) Sh — 0 
Vd =~ rel 
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60 


65 
EO 3 


70 
104 
80 
FO> 
90 


106 


20 
og 


eat. Gir 20) GOUTO «65 
TBAR = Q.DO 
DO 60 KK = IK,IKE 


IF (IX(KK,J) .EQ. 0) TBAR = TBAR + 1.D0 
0) TBAR = 


me (f.Gt. 1.D—<4 .AND. IX(KK,dJ) .NE. 
pee sey KK, J ) 

CONTINUE 

IF(TBAR .LE. 1.D-4) GO TO 70 

TBAR = TBAR/ANN 

CALL ANEWT (TBAR, PAST, XBAR,T,PM2) 
HATM1 = XBAR/PM2 


IF (DABS (HATM2 - PM2) .LE. 1.D-4) GO TO 80 


HATM2 = PM2 

Cee ro 30) 

WRITE (6,103) AM1,AM2,BM1,BM2,TI 
HORMAT('O',114,4F3.3,T49,F5.4,T56, ' 
BO N ye ly ae) 

Ger To 90 


himetoco, 104) AMI, AM2,BM),BM2,12,0,8M1,BM2 ,JJJ 


TBAR + 


Meri Oo La Aro .3 5049, bora oo er .4,2F8.3,15, ° %") 


Gee TO 90 


Mimi (6,105) AML AM2,SMl,BM2Z,TL,1,HAIM] ,PM2,J77 
Bemtinn('O ,T1l4,4F3.3,1749 ,65.4,1756,F5.4,2F38.3,215) 


CONTINUE 
WRITE (6,106) 
FORMAT ('1') 
Sor 

END 


SUBROUTINE ANEWT (TBAR,YIN, XBAR, 


IMPLICIT REAL*8 (A-H,0-Z) 
COMMON HATM1,HATM2 

TL = DLOG(TBAR) 

DO 20 KKK = 1,20 

GY = DEXP (YIN* (T-1.D0) ) 

FY = YIN*TL + XBAR*(1.D0 - GY) 
FPY = TL - XBAR*GY*(T - 1.D0) 
YOUT = YIN - FY/FPY 

IF (DABS(YOUT - YIN) .LE. 1.D-4) GO 
YIN = YOUT 

CONTINUE 

RETURN 

END 


SUBROUTINE ODT. (7) 


BMP CitenbAL~ 3 (A-H,O=Z) 
COMMON HATM1,HATM2 


FLAG = 0 
THETA s=—9299 .b-—4 
EZ = 30.80 


FIND A COARSE BRACKET AROUND THETA: 


70 


fe GO E) 


TO. 99 


(A,B) 





Js 


Z0 


ja, 


36 


40 


THETA = THETA - 1.D-2 
Pra p ia. Lt .0.D0) GO-TO 15 
i = BPE (THETA) 

Patt ae. EZ GO TO 1s 


nee Bo 
GO TO 5 
A = THETA 


Pao. oT. 62. bD0) A =" O2nDe 

B = THETA + 2.D-2 

IF(B .GT. 9999.D-4) B = 9999.D-4 

R = (-1.D0 + DSQORT(5.D0))/2.D0 

Al = (1.D0 - R)*(B - A) +A 

Bl = R*(B - A) +A 

EAL = EFF (AL) 

EB] = EFF(B1) 

IF (DABS (EAl-EBl1) .LT. 1.D-15) GO TO 30 
IF (EAL .GT. EB1) B= 

DEe4{EAL .LT. EB1) A = 
coero 20 

me (Al + Bl) /2.D0 
CHECK FOR MAX FROM THE LEFT END 
MeCELAG .EO. 1) GO TO 40 

IF(T .LT. 1.D-1) RETURN 


Bl 
Al 


TSTAR = T 
ESTAR = EBL 
THETA = 0.DOQ 
E2 = 0.D0 


BETA = THETA + 1.D-Z 
Peta oeGl. 99.b-2) RETURN 
ET = EFF (THETA) 

fie el. £2) GO TO 36 
B2 = ET 

GO TO 35 

eee TAEDA = 2.D=2 

mA Gl.) 0.D0) A = 0.D0 
B = THETA 

FLAG = lL 

co TO 20 

Bei porak «.LT?.EB1L) RETURN 
T = TSTAR 

RETURN 

END 


FUNCTION EFF (Q) 


IMPLICIT REAL*8 (A-H,0-Z) 
COMMON HATM1,HATM2 

AM1 = HATML 

AM2 = HATM2 

DEXP (AM2*Q - AM2) 

Bec AME~ DOs). 
F2 = DEXP(AM2*Q*Q - AM2) 

CR = Wiay.ce (elie (CIE BIO) 3 182 9) 


seh 
G 


Ge 





EN = 1220) (AMZ- Gan ee DO) 

FN = AM1L*G*G*FN*FN 

Oe Cara 2AM nM 60@~ bn 0 1 (O7r es 0) 
FD = AM2*((1.D0 + AM2)¥*G2 - G*G*FD) 

EFF = FN/FD 

RETURN 

END 


eZ 





COMPUTER PROGRAM 4 


C RANDOM NUMBER GENERATION FROM NEYMAN'S TYPE A 
C DISTRIBUTION 
DIMENSION AD(50) 
DATA ISEED1,ISEED2,ISORT,KSIZE/431219,187648,0,50/ 
CALL OVE LOW 
DOwLO ML = 1,9 
PM1 = FLOAT(M1) 
CALL LPOIS1(ISEED1,AD,KSIZE,PM1,ISORT) 
WRITE (6,1) 
iL FORMAT ('1',///////) 
WRITE (6,2) 
2 FORMAT (' ',14X,'RANDOM NOS. FROM THE NEYMAN TYPE A 
* DISTRIBUTION’ ) 
WRITE (6,3) 
3 FORMAT ('0O' ,//,12X'Ml1 VALUE' ,T40,'M2 VALUES') 
CALL DGHTER(AD,ISEED2,M1,KSIZE) 
10 CONTINUE 
STOP 
END 


SUBROUTINE DGHTER(AD, 1SEED2,M1) 
DUMENSTON AD(50),Y( 100), IX(30, 3) 
Wai o,99) (0 ,J—1,9) 

29 PORMAT(*O',18X,915,//) 
DO 200 K = 1,50 
ISIZE = INT(AD(K) ) 
Meat Size «NE. 0) GO TO 53 


BO 50 M2 = 1,9 
50 IX(K,M2) = 0 
GO TO 55 


a6: po L100 M2°>= 1,9 
PM2 = FLOAT (M2) 
See beOlsl (ISEED2Z,Y,usige, ez, 0) 
IX (K,M2) Q 
Bo LOO 1 lone 
100 eck, M2)o= IX(K,M2) = INE CYT) 
a5 Nieetee( Go, Ol) Mig (ix (hee) 
ro 1 BORMAE Ge @ pLOX i> ,3X%7,515) 
200 CONTINUE 
RETURN 
END 


; ou di 


ya 





APPENDIX B 


RANDOM NUMBERS FROM THE 
NEYMAN TYPE A DISTRIBUTION 


Computer program 4 was used to generate the eighty-one 
sets of random numbers for the numerical examples of Table 3. 
The parameter range used was my and m. from 1 to 9 in incre 
ments of 1 and the sample size was 50 in all cases. The 


following pages contain the data sets generated. 
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RANDOM NOS. FROM THE NEYI AN TYPE A DISTRIBUT ICN 


M2 VALUES 


M1 VALUE 


GO TPOQVDOAD DOGO PBONGOF OM NATFOONANOR OMONNOANODOON + 00M 
ad Sri “s~ ANH AICS ex do -4 


MAINO OSORONOLOOIPHRIOOAN OF SOMBONN DOWD TDODQOVINEONOV4.00IN— 
s+ st = mt ed mt tt 4 


EL OD OOD I ODI SOY OSMOSIS MMOMAN-ONNONYVOFOOMOO in 
ool se eee ho oe Po oad — 


FOSS GOOF OMODOMOOVOT 4 DONG VOANINODMONTNITOS.QOOr meomnin 
eS Oe me) NOR ORO GOONIES N 
IL DODAGL IMAI OL ONG OGRA GOIN WA AMAINIM Din FODOOMO DOIN 
te SOON AAO AMO CGH CG aN ter 


MAMNGOTOVOMOVONAOMOOS ANHOFOO OM HOWMNT HD NODNO QOOdR EMA 


AOOMODVIOOVIODOVOOS NN GAO MM AN AOIOOONFHOFHONOMTOOFHOYOOQOMNG0N 


UN eh dF ng ak md pad ed Ty Nm ad ad ng a a ed dd Gd eh ay de hy a ad nd Og et ed md 


72 





~~ aa 


FROM THE NEYMAN TYPE A OISTRIBUTION 


NCS. 


RANDOM 


M2 VALUES 


Ml VALUE 


ee Oe nN ee ON OO OD OF NOMP ONS ONTO 
NWN aeaNqaqmf MM ANSE NOt Mt ANNA | ON Na oA t+ 


DO TONONYNDHOMAOSFOROMONMAT RIND D+O0MOF-OMOKrVONONTOODOF OO 
N MQ ett SteNNA NN ANS SB Qt OS atta A OMe Neti A 


DOVOSO QODANMDOOMOAT OONS FADS N SF ODSSFMNOANIDAATATFNBWDANDAAS ODDO 
NO AANAANN A ANG N mt Net NANA a ound NA a — 


OOO Oe Se ne tig cna a COO NOON TS DO MAO DIA OINAS 
rw SALaedtNay et aN, NA NS eet Net 


ears | See a ne ee SON PON OM OTOP NS 
NN atin ate etq toe NA teed ct 


Se Se Oe oa ee OO a ea OOS Ae OM OOM anoms 2 
mse NN tt nt mir tet 


eee Aah) eg On ce oe ha OS INOS NO SS 
a4 4 aN] onal 


OOWMOWO MOR LALIVRE OHO MnyOOrn]A DDODOANOMAR MND ALS F VOOM NOMAD 
= 


VOFOTOMIOONVNAOQNINOVNANDVATDD-F OMAN ANN DOONOINNNQOQANNAAOQOrO 


PNENLAIIEN AIO OTIEN AEN TAIT EN AION AT CIO EN CUENTA OU IEN QI EN TATOO QUOI AI OI ATA EN AION AIO ON IOI ONOT 
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FROM THE NEYMAN TYPE A DISTRIBUTICN 


RANDOM NOS. 


M2 VALUES 


Ml VALUE 


et OOO ON Co UN 2 Orel 2 COUN ne Gye coal i NT 0 DO ODO OAT ODOM 
MrNaoanneo FHKAMANNEN GQ NOAM NRenNoNt At toereae me 


SO rene ee Oe oem 
OT AAD aA MO SAUL EAM GANNAIN HRN O A HANAN A QyamNmOee NO 


DIANL HSNO, FOOD DADIN™M ANMNAqHOADDDAN OOFINDHRANDOFA+90 OMA FoF O 
NMA Faas MO SANA AAOAN CQ ties QC et HNN aN Jas 


ee Oe eo ee DADQHAD HOCOMAM DMN QADOODOAHP-&—ODOrO 
CO SF Ag Rt et 4 Teqrtayqr*tqugdeteq ama aN et NaN aetqQQyge ACN 


NP DAM DP SFAIONOLIN BINS AO MMDOATUOAM OVD VE NNODOAOODONAD QNONOMNM 
NN edeHtN seq, NS MQrtaqgeetqNatae —-a cH mt ont gg NOY NAN OR NO 


OED OOO ANG ag Oro OR, ee ee Na aon 
4aN4 WC SA tte IN rt mt et nd mat HAIN YY 


Ne Nae re IN ae COR OOS ST VOI OS IMO OF 
Aedes ed — nad Ait et Ate _ m4 mt NI 


SSSR eset) ee eg GO A ot ee 
etree 


PHET OSS MMOMOADM HMMNADTANAQVHNG DO NNMGANN SE OONOMFUONE VOM 


OY YTV LEY EY EVV CY VDE} EY EVH EY CH. OEY.) CVA ET HEY EY EV. OY COD EF FE EY. OD. EYED EDF V EYE FVD EYED FY CTY FD ECVE CTY ETT OD 
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TYPE A DISTRIBUT IDN 


FROM THE NEYMAN 


RANDOM NOS. 


M2 VALUES 


ML VALUE 


FODIMOOMO AAT ANODAD QHMAOIAINOM OF AIAINO DINIF INDO + FOODAON— 
FOO SOE OR NOP NIN FLAN AANO OK mt tH EOI ANION St NIL ALO et St LEN CO OI 


ROWAN AIA HOO OD FAD AY SN O DOP TODOS DAM AT CTP AP PY AP YE OPM OT O 
AYN SE tC PF ON AIAN EN OPO Oe EAE MAOANN SNA NET OE OTN 


ILM IR OID KI MDA NOM $F ON ADO MODONMM AED DOK OMAING AIO OOM 0 $0 + 
NaN AIO Maa Noah Naa WON HQ staqQnian Nanas atiag NAILAN AIS 


FOO MD MOAR MOAODN HODNAD LO VATE FOO AOR AN AONNO ND OR OM QOrtes 
TaN tsa MNT MOM ANNA at ANNA et tat HAIN TAM TY CHLAN 4 


MMO ONG IO TOMNFOHODMOOTDAEROVNANMMS VF NOODONDT NNO OO DHUE OD 
aNdaaNcag neta Qetet NaN aaNet tile yeasts aNest ate Ost MAM 


ARNO ROR NNO ATR OT NAD MOPFOMIAALOIE OM Qt ANDDO0POAN AMM MOM O+ 
QING Neamt ANMNNN ATMs AM Nest acetic aQyateam) AN ati 


actA, aReaN aa apaNet tay at at rt HA oa a ae aN ne ha eee) 


DONYODOEMINMVNGONDON S DHOBRNAOT ADO OM ALD OMIM FABIIR AT AAS ON 
a = — HA ete mj CJ —_ -4 not mt end 


— ee | 


FPP PTH TP HT LST PHS TH TH TF HT HS LMS SHSM HST HAMAS HHH TAIL HT 
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NEYMAN TYPE A OISTRIBUT ION 


ets 


FROM 


RANDOM NOS. 


M2 VALUES 


M1 VALUE 


DONO ODDDOFMOROQD OR TR Ot NN Mh unl MAMA ME AGH OF DOOR OR tFINNOCIM,.O 
MminOodoina $Onyo 4 FOTRMAMNMOA ONS RS ONMADOF MOF OINME QATR HOHE oO 


WQLATIO SOLO CIC UN DO FONAQO HM AAE EH ON DOP DADWINDMA MQWOO FUAk- Pink OMTUE 
TIANYVINS Fay DRAIN SB MOT R AMS TPO OOM AM ST VOTO NH MOADHNME 
| 


D amaAD OM AO OM Ff WIN AY. FO AN NIA NI 9 119 0 ND ON PRIN DLA OOD SLA ANI ah Doo 
OAS Pean FOV P 2 OFAN OPIS FOF FAIA DOIN MNS OUIALINAL MNnADo ato 


DOT ONAOKR..DQDHNOGDNUSTAMOENODDAMN SQ reaNn st ODOR MOHIND ODANiNn 


Ot oro 
mam Fen Net MoamM NOISES HF NANFAN ONION FENATMMMOAN NANA AMY 


DRAAMQODAHOWMOADA FOR DM AORN D AOD QIN Feattinat Om VUNGHODOVONMO OMEN 
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NIP DEDOAMQOOMUOW FAROE ES ODOR. DOE FAIA FRO OMOAMONODO Far nintoon sto 
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DINA ROME ON AN OQ QyARQIN MO MN OVI NONE © SOR FH OA AIO DOF LAM OFFER HOO 
NamaNal Nay Sa Stee AQMAAaA NRA SIRO HON ON RS OHM 


SO 2 Oe Oe eee Sie Ten OO OO NOS Sar 
aot ay aN | mt =| —4 a ae AAs Set ON 


FS SI Cl a Su ar Ae 0 ame ine SN OND Oot 
— 


LOY LED EID) EA EG UD ED Uy ETUDE VED 9 8 ny UV EDT ER UD EVE) COTE ER ED COVE EDV EG LD SER UNV EV ED EQ EDV iy LAL UY 


ie, 





TYCE feo oe Een 


FROM THE NEYYAN 


RANDCGM NCS. 


W2 Vv AEGES 


ML VALUE 
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FROM THE NEY4AN TYPE A DISTRIBUTION 


RANDOM NCS. 


M2 VAIWES 


ML VALUE 
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RANDOM NCS. FROM THE NEYMAN TYPE A DISTRIBUTION 


M2 VALUES 


ML VALUE 
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FROM THE NEYMAN TYPE A DISTRIBUT ICN 


RANDOM NOS. 


M2 VALUES 


ML VALUE 
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